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In this paper we introduce a generalization to the algebraic Bender-Wu recursion relation for the 
eigenvalues and the eigenfunctions of the anharmonic oscillator. We extend this well known formal- 
ism to the time-dependent quantum statistical Schrodinger equation, thus obtaining the imaginary- 
time evolution amplitude by solving a recursive set of ordinary differential equations. This approach 
enables us to evaluate global and local quantum statistical quantities of the anharmonic oscillator 
to much higher orders than by evaluating Feynman diagrams. We probe our perturbative results by 
deriving a perturbative expression for the free energy which is then subject to variational perturba- 
tion theory as developed by Kleinert, yielding convergent results for the free energy for all values of 
the coupling strength. 

PACS numbers: 05.30.-d 



I. INTRODUCTION 



Most physical problems can only be solved by approx- 
imation methods. One of them is perturbation theory 
which yields weak-coupling expansions. Unfortunately, 
they often do not converge. 

The ground state energy of the anharmonic oscillator 
is the simplest example where this phenomenon can 
be studied. Algebraic recursion relations a la Bender 
and Wu jlj] yield perturbation series for the eigenvalues 
(energy) and eigenfunctions (wave functions) of the 
time-independent Schrodinger equation up to arbitrarily 
high orders. In Ref. ^ the calculation was extended to 
250th order. The Bender-Wu recursion relation yields 
a power series for the anharmonic part of the wave 
function both in the coupling strength g and in the coor- 
dinate x. The power series in x can be cut off naturally 
by comparing the recursive results with those obtained 
from evaluating Feynman diagrams. The resulting 
weak-coupling series for the ground state energy does 
not converge for any value of the coupling strength. This 
paper deals with both problems: Obtaining high-order 
perturbation expressions and making them converge for 
all values of the coupling strength. It is organized as 
follows: 

In Sec. II we perturbatively evaluate the path inte- 
gral representation for the imaginary-time evolution 
amplitude of the anharmonic oscillator by means of 
a generalized Wick's theorem p, [ 



In Sec. Ill 



we 

represent the first-order results diagrammatically. Doing 
so, we demonstrate that the algebraic computational 
cost is very high for the diagrammatic approach. We 
also obtain a cross check for the results which are 
derived from a differential recursion relation in the sub- 



computational cost we introduce a strategy to exploit 
the symmetry property of the imaginary-time evolution 
amplitude in Sec. [v| thus laying the foundation for 
our high-order results. In Sec. VI we combine the 



resulting algebraic recursion relation with the original 
differential recursion relation, thus generalizing the 
Bender-Wu approach From our perturbative results 
for the imaginary-time evolution amplitude we then 
gain a perturbation expression fo r th e free energy of 
the anharmonic oscillator in S ec. |VII| which we check 
again diagrammatically in Sec. |Vinj Th e perturbative 

by means of 
intermediate 



IX 



sequent Sec. [V. In order to cut down on the algebraic 



results are then re-summed in Sec 
variational perturbation theory || for 
coupling .9 = 1 for which the usual weak-coupling series 
would diverge. This theory is a systematic extension 
of a simple variational approach, first developed by 
Feynman and Kleinert in the path integral formalism. 
Feynman introduced the path integral formalism as a 
quantization regulation, that represents the operator 
properties of quantum physics by fluctuations of the 
dynamical variables ^ ?j. By extending analytically 
real time to imaginary time, also quantum statistical 
quantities can be obtained by summing over quantum 
mechanical and thermal fluctuations with the help of 
path integrals [pi H|. In order to evaluate the path 
integral for the free energy approximatively, Feynman 
and Kleinert developed a variational method in 1986 ||. 
It replaces the relevant system by the exactly solvable 
harmonic oscillator whose frequency becomes a varia- 
tional parameter which has to be optimized. Starting 
with Ref. |l(J, this method has been systematically 
extended by Kleinert to higher orders || [H]]. It is now 
known as variational perturbation theory and yields 
results for all temperatures and all coupling strengths. 
In Sec. [x] we extend this procedure to higher orders of 
the free energy and cross check the results in Sec. XI, 
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II. PATH INTEGRAL REPRESENTATION 

The path integral representation for the imaginary- 
time evolution amplitude of a particle of mass M moving 
in a one dimensional potential V(x) reads || 



(xt hf}\x a 0) 



x(hf3)—Xb 



x(0)=x a 

M 



hp 



(h 



x 2 (t) + V(x(t)) 



For the anharmonic oscillator potential 



V{x) 



M 



cu 2 x 2 + gx A 



(1) 



(2) 



the imaginary-time evolution amplitude (|]J) can be ex- 
panded in powers of the coupling constant g. Thus we 
obtain the perturbation series 



{x b h/3\x a 0) = (x b hj3\x a 0) u 



(3) 



fit) 



dTl(x A (Tl)) u 



where we have introduced the harmonic imaginary-time 
evolution amplitude 



(xb hf}\x a 0) u 



x exp ■ 



x(hl3)—xb 



I f HP 

dr 



ft Jo 



Vx 



M .,, . M , ' 



(4) 



and the harmonic path expectation value for an arbitrary 
functional F\x\: 



(F[x]) u 



1 



(x b hfi\x a 0)o, J x{0)=x 



x(h(3)=x b 



VxFlx] 



(5) 



x exp 



Tr3 



dr 



M .,, . M 2 2/ s 

Y^ 2 M + yw 2 s 2 (r) 



The latter is evaluated with the help of the generating 
functional for the harmonic oscillator, whose path inte- 
gral representation reads 



(x b hf3\x a 0) 



x{hj3)—X}j 



2?ccexp 



s(0)= 



1 

ft Jo 



dr 



M -2, s M 2 2, ^ ■ f \ i \ 

Y x ( T ) + y w x ( T ) -J\ T ) X \ T ) 



leading to |j| 

(x 6 ft/?|a; a 0)a,|j] = (x b h/3\x a 0) a 



(6) 



(7) 



x exp 
1 



h8 



dri x c \(Tijj{Ti, 



ii 



+ 



h/3 ph0 

dn dr 2 G( D '(r 1 ,r2)j(ri)j(T 2 ) 
Jo 



with the harmonic imaginary-time evolution amplitude 

(x b h(3\x a 0) 



Mlo 



exp 



27rft sinh h(3uo 

-[{x a + x 2 ) cosh hf3ui — 2x a Xb] 



(8) 



2ft sinh ft/3u; 

In equation (Q) we have introduced the classical path 

. . x a sinh(ft/3 — r)u> + Xb sinhwr 
* cl(T) S shrhft^ ' (9) 

and the Dirichlet Green's function 
ft 1 



G( D )(r 1; r 2 ) 



(10) 



Mui sinh ft/3w 

— t 2 ) sinh(ft/3 — ti)w sinh WT2 
+ 0(t"2 — T i) sinh(ft/3 — t 2 )ciJ sinh wti] . 



We follow Ref. j| |[ and evaluate harmonic path expec- 
tation values of polynomials in x arising from the gener- 
ating functional (0) according to Wick's theorem. Let us 
illustrate the procedure to reduce the power of polyno- 
mials by the example of the harmonic path expectation 
value {x n {T l )x m {T 2 )) u} : 

(i) Contracting x{t\) with £™ -1 (ti) and x m (r 2 ) leads 
to Green's functions G^^n) and G (d) (ti,t 2 ) 
with multiplicity n — 1 and m, respectively. 
The rest of the polynomial remains within the 
harmonic path expectation value, leading to 
(x n - 2 (n)x m (T 2 )) u and (i^-^n)!™- 1 ^))^ 

(ii) If n > 1, extract one x(ti) from the path 
expectation value giving x c \(ti) multiplied by 

{x n - x {rx)x m {j2))o J . 

(hi) Add the terms from (i) and (ii). 

(iv) Repeat the previous steps until only products of 
path expectation values (x(ti))^ = a; c i(ri) remain. 

With the help of this procedure, we obtain to first order 

(z 4 (ti)) w = x^(n) 
+ 6 4( n ) (ti, ti) + 3 G^ 2 ( n , n) . (n) 



III. FEYNMAN DIAGRAMS 

These contractions can be illustrated by Feynman di- 
agrams with the following rules: A vertex represents the 
integration over t 



dr . 



(12) 
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a line denotes the Dirichlet Green's function 



and the second diagram reduces to 



(13) 



and a cross or a "current" pictures a classical path 



Zcl(n) • 



(14) 



Inserting the harmonic path expectation value ( |ll| ) into 
the perturbation expansion (|§|) leads in first order to the 
diagrams 



hi) 



dTi(a; 4 (Ti)) u 



(15) 



;+ 6 x ^ x + 3 



We now evaluate the first-order Feynman diagrams in 
( |l5| ) for finite temperatures and arbitrary x a ,Xb- Thus 
we will get a first-order result for the imaginary-time evo- 
lution amplitude in (pi). The first diagram leads to 



1 



[(x^ + xl) (sinh4/i/3w 



(16) 



32Mw 2 sinh 3 hfiu 
+ 9 sinh hfiuj — 12h(3uj cosh H(3oj) 
+ x a Xb (— 12 smh2hf3LU + I6h[3uj 
+ Eh/3u cosh 2h(3bj)] , 



[(xl + xl) (sinh3/i/3w 



(17) 



whereas the last diagram turns out to be 



' 16M 2 u 3 sinh 2 Hf3uj 
■ihfiLu + 2hf3Lu cosh2h/3u)) . 



-3 sinh 2h(3oj 



(18) 



32a; sinh h/3u> 

- 8sinh2fi,/3w + 12H/3uj) 

+ (x^Xb + x a Xb) (4 sinh 2>hf3u) + 36 sinh h(3uo 

- 48h(3uj cosh H(3w) So all in all we get tne following first-order result for the 
+ x\x\ (—36 sinh 2h(3uj + A8h(3uj + 24hf3ui cosh 2hj3uj)\ , imaginary-time evolution amplitude 



(x b h(3\x a 0) = (x b h[3\x a 0) u 
h 2 



9 



x 1-i 



h { M 2 iu 3 sinh 2 h/3u> 

h 



Muj 2 sinh 3 H/3u) 
9 



9 3 3 

sinh 2h/3uj H — ft/3o; H — fi,/3o; cosh 2fi,/3w 

16 4 8 

27 



(x a + x b ) — sinhihf3u> + — smkhflu — — h[3ui cosh Hf3u> 



16 



16 



.'',,.17, j — — sinh2h[5u> + 3h/3uj + —h(3Lu cosh2/i/3u; 



1 



to sinh hf3u 



(x a + xl) ( — sinh4S/3w - - sinh2?i/3w + -h/3w 



+ (x 3 Xb + x a xi) — sinhZhf3u} H — smhhf3uj h(3uj cosh h(3uj 

\ 8 8 2 

/ 9 3 3 

+ x 2 xl I — — sinh2/j/3w + —Hflu + —H/3u; cosh 2hf3uj 
V 8 2 4 



(19) 



r 



The imaginary-time evolution amplitude thus has the 
time reversal behaviour 



while it is known that the imaginary-time evolution 
amplitude is real for one-dimensional problems. 



(xb HP\x a 0) = (x a h(3\xb 0)* , 



(20) 
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IV. PARTIAL DIFFERENTIAL EQUATION 

Consider the Schrodinger equation for the real-time 
evolution amplitude 



d 

lh—(x b t\x a 0) 

at 

h 2 d 2 

~ ~ 2M dx 2 



(21) 



(Xb t\x a 0) + V(x b ) (x b t\x a 0) . 



In order to get a corresponding quantum statistical 
Schrodinger equation we now have to change from real 
time to imaginary time, i.e. we have to perform the Wick 
rotation t — ► —ir. Thus the Schrodinger equation (pi 
becomes 



d 

h—(x b T\x a o) 

OT 



For both the real and the imaginary-time evolution am- 
plitude the initial condition reads 



(xb 0\x a 0) = 8(xb - x a ) 



(23) 



Plugging the anharmonic oscillator potential (g) into the 
Schrodinger equation (p3) we finally get 



8 



h 2 d 2 



M 



dr 2M dx 2 



2 2 4 
^ %b ~ 9 X b 



(24) 



(x b r\x a 0) = . 
Making the ansatz 

(x b r\x a 0) = (x b r\x a 0) w A(x b , x a , r) , (25) 



where (xb r\x a 0)^ is the harmonic imaginary-time evo- 
lution amplitude (^), we conclude from ( p4|) a partial 
differential equation for A(xb, x a , r): 



d_ 

Ot 



h d 2 

2Mdxi 



Xb cosh lot — x a d 



sinhn 



dxb h b 



A(xb,x a ,r) = 0. 



(26) 



We now choose our ansatz for A(xb,x a ,r) by introduc- 
ing three expansions in g, in x a and in Xb, respectively. 
Also we take out the factor sinh - ' lot, such that the ordi- 
nary differential equations for the expansion coefficients 
become as simple as possible: 



oo 2n 2k 



A(x b ,x a ,T) = J2J2Y,9 n - 



» 

2k\l 



(r) 



,2k-ll 



n=0 k=0 1=0 



sinh lot 



(27) 



In order to obtain the unperturbed result A(xb, x a ,r) = 1 



+ V(x b )(x b T\x a 0). (22) for g 



we need c^ (t) 



1. The superscript n in 
equation ( |27| ) denotes the perturbative order, whereas 2k 
counts the (even) powers of the various products x l a x\. 
The summations over the coordinates be trun- 

cated at k — An, because we learn from Feynman di- 
agrammatic considerations that the diagram with the 
most currents x in the nth order looks like 



(28) 





Inserting the new ansatz (27) into the Schrodinger equa- 
tion ( p6] ) and arranging the indices in such a way that 
each term is proportional to x 2k ~ l x l b , we get for the dif- 
ferent powers of g and for n > 0: 



J 



EE Ss£^ - £ E E" + m + d^SS^ 

fc=0 /— /c— — 1 (=— 2 



2n 2fe-l J n ) (-\ 1 2n 2fe+4 J^-l) / \ 

fc=(H=-l Smh UT h k=2 1=4 Smh UT 



Thus the sums over fc and over I collapse and we deter- which is solved by 

» 

J 2k\l 



mine the master equation for our coefficients cl?j , (t) 



t (n) M c " fc ; (r) = (l + 2)(l + l)— dr ^ - (31) 
9c 2fc|i( T ) _ , „ - fi c 2fc+2|;+2( r ) ' 2M ^ smh wr 

9r -^ + 2K/ + iJ 2Af S inh 2 ^r /• 42l +1 W 

(„) + (/ + l)u / 

c ?i.i/j_A T J 1 ii „ J smh wt 

smh wr ft - W ^4fc-iUw sinh WT + d S 
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Here the d^j; denote the integration constants which are 
fixed by applying the initial condition 



(ii) Drop all terms containing a ci?], (r) with Z > 2fc. 



lim 

T-.0 



c (n ? (V) 



sinh' lot 



< oo . 



(32) 



(iii) Neglect all terms containing a c^j ; (t) with any neg- 
ative indices k and Z. 



However the above master equation ( p0[ ) is not valid for 
all fc and Z. Therefore we now introduce a set of empirical 
rules telling us which of the coefficients c^j ; (r) have to 
be dropped once we write down (^lj) for any order n: 

(i) Drop all terms containing a c^Ju^r) where 2/c > 4rt. 



To convince the reader that equation ( |3l| ) together with 
this procedure leads to the correct results we now reob- 
tain our first-order result from (|l9|). To that end we set 
n = l, such that k runs from to 2 and Z from to 4. 
Fixing k = 2 and counting down from I = 4 to Z = we 
get 



c (1) (t 

C 4|4V T 

c (1) fr 

C 4 3V T 



c (1) (t 

C 4|2V 



c (1) fr 



c (1) fr 



- y drc^ (t) sinh 4 lot + d 



ft 

4lo I dr 



3lo / dr 



sinh ut 

sinh 2 ut 

cjg(r) 
sinh 2 ut 

sinh cjT 

Correspondingly for fc = 1 we obtain 
6ft 



= 2lo dr 



lo dr 



d W = 

"4 3 



d W = 
"4 1 



ijj = — [ — sinh4tjr sinh2o;T H — sinhwr 

4 I 4 ftw 132 4 8 



1 



1 9 3 

.— sinh Slot +— sinh lot — — cosh ut I , 
ftw smh cjt \ 8 8 2 

1 / 9 . , n 3 3 

2 — — smh 2lot + -wt + —lot cosh 2wr ] , 

fiwsinh wt V 8 2 4 



1 (1 . , o 9 . , 3 

: — g — smh3tJT H — sinhwr wrcoshwr 

ftu; sinh lot V 8 8 2 



1 



' huj sinh wt V 32 



113 
— sinh4tJT — — sinh2ojr + — sinhwr ) . 
4 8 



c«(r) 



ill 



dr 



48 w 



4Sw 



sinh 2 lot 

(1)/ \ 

M j sinh wt 



3 27 , 9 

,i, — ,,„.., . — smh3cjrH — - smhur — -cjrcoshwr 

2 I 2 Afw 2 sinhwr I 16 16 4 



2w / dr 



sinh 2 (r) 



7(1) 



M lo 2 sinh 2 lot 



9 3 
- — sinh 2lot + 3ojt H — lot cosh 2cjt 
4 2 



1 



3 



ft T c[]1(t) r 

— / ctr — hw / dr — i-s — 

M J sinh 2 wr J sinh 2 (r) 



l 2\0 



27 . 9 

■ -sinh3u;r-| sinhwr ojrcoshojr 

Hu> 2 smh 3 lot V 16 16 4 



(33) 
(34) 
(35) 
(36) 
(37) 

(38) 



(39) 



(40) 



Finally for k = we get the equation 



ft / , l -2|2 

M ■ dT 



c«(r) 



sinh wr 



d (1) - 



9 • , „ 3 3 

. smh2cjrH — lot H — a>rcosh2ajr 

M 2 w 3 sinh 2 wr V 16 4 



(41) 



r 



The path of recursion which follows from this procedure V. EXPLOITING THE SYMMETRIES 

is shown in Fig. |l|. 

As seen above we already have to solve nine ordinary 
differential equations for the first-order imaginary-time 



6 



evolution amplitude, 
integrals to solve is 



For any order n the number p of 



(|43|). Comparing equation ( pCf ) for k = 2, 1 = 4 and 
k = 2,1 = we then obtain an algebraic equation for 



c^(t). The knowledge of c^(r) gives us c^j^r) be- 



„(i) 



2n+l 

^ (2j - 1) = 4n 2 + 4n + 1 . 



(42) 



» 



(r) 



(n) 
C 2k\2k~l 



cause of the symmetry ( ]43| ) and by comparing (30) this 
time for k = 2, 1 = 3 on the one hand and k = 2, 1 = 1 
on the other hand we are left with an algebraic equation 
Due to the time reversal behaviour (^0|) , the coefficients for c£2 (r) . Thus we get all the coefficients for k = 2 only 
C 2fc|/( T ) snow a symmetry, namely: by solving one differential equation, namely the one for 

c£m(t). For k = 1 the procedure is similar, k = only 

generates one coefficient anyway, namely Cq\ (t), which 
still has to be solved by evaluating one integral. The new 
path of recursion is shown in Fig. |[ 
So finally three out of the nine first-order coefficients are 
obtained by integration, three more are clear for symme- 
try reasons and three come from an algebraic recursion. 
We now generalize the algebraic part of our recursion. 
Consider again the symmetry property ( f43| ) . Differentia- 
tion on both sides yields 



sinh' lot 



sinh 2 ' £ 1 lot 



(43) 



Exploiting the symmetry ( (43[ ) , we can cut down the num- 
ber (ff|) considerably. At first sight it is reduced to 



2n+l 



p' = J = 2r 



3n + 1 . 



(44) 



so there are only six integrals left for the first order. 
But we can go even further. Employing these symme- 
tries we can eventually change almost all recursive dif- 
ferential equations into purely algebraic ones leaving only 
p" = (2n + 1) integrations. So for the first order we are 
left with three integrations only, namely with equations 



2k\l 



(r) 



sinh 



2k-2l 



@8D, and (|4lD 



These coefficients c^{|(r), cij^r) 



-2(fc — V)oj coshwr 



LOT 

» 

J 2fc|2/c-i 



sinh 



2fe-2;+i 



(45) 



4|4V <" "-2|2 V 

and Cpjp(r) are integrated recursively. The other coeffi- 
cients can then be obtained algebraically: Once we have 
cju^r) we also know c^fY) because of the symmetry 



Now we substitute for the two partial derivatives accord- 
ing to equation (|30|). Solving for the (I + l)-st coefficient 
and shifting the index I down by one we obtain 



J 



(n-l) / \ 

(„) _ {I + 1)H (») C 2k _ i{l _ 5 (T) 6 

C 2k\M> - 2Mu) C 2k+2\l + A T )+ ^ SUIH WT 

(*V) / \ (n) / \ 

, (2k-l + 3)(2k-l + 2)hc 2k+2l2k _ l+3 (T) 2k -I + 2 c 2k\2k-i+2( T > 



2MloI sinh 2fe - 2i+2 lot I sinh 2k - 2l+2 lot 

1 c 2k-4\2k-i-3y T ) {2k -21 + 2) cosh lot c 2k\2k-i+\ y T ) 



fuol sinh 2fe - 2 '- 4 W T J sinh 2fe - 2i+1 



which is the algebraic recursion relation for any non- 
diagonal coefficient <4fc|i( r ) with < I < k. (The coef- 
ficients with k < I < 2k are then clear for symmetry 
reasons.) The diagonal coefficients <4fej 2 jfe( T ) have to 
be integrated. 



evaluated by integrating the differential recursive equa- 
tion, we can even further simplify the solution ( |3l| ) to our 
master equation (|3^). We only need it for the diagonal 
coefficients, for which Z + l = 2k + 1 is always greater than 
2k. And according to our index rule (ii), coefficients of 
the shape c^, 2k+1 have to be neglected. We get 



VI. COMBINED DIFFERENTIAL AND 



ALGEBRAIC EQUATION % f 4™ > , Jr) 

4fe2fcW = ( 2k + 2 )( 2fc + 1 )wT7 rfr 2fc+2 l 2fc+2V 



2M J sinh lot 

We now combine the differential recursion with the al- 1 f (n-l) ,4 ,(n) 

gebraic one. As only the diagonal coefficients have to be hi 2k-4\2k-4\ T ) smn UT + "2fe|2fe • ( 47 ) 
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01234 / 01234 I 



FIG. 1: This diagram depicts the path of recursion for n = 1. 
We start in the top right hand side corner, which is to be 
identified with the coefficient cffi and follow the arrows until 
reaching the bottom left hand side corner with the coefficient 

c (1) 



Index rules (i) and (iii) still have to be applied, k runs 
from to 2n. 

Let us quickly summarize the combined differential and 
algebraic recursion relation considering the first order as 
an example. Fig. |^ shows all first-order coefficients for 
the imaginary time evolution amplitude. Each coefficient 
is represented by a little circle. Now the coefficients on 
the diagonal line 2k — I have to be obtained by referring 
to equation ( {47] ) together with rules (i) and (iii). These 
two rules tell us which of the coefficients either from the 
the same order n or from the previous order n — 1 have 
to be integrated and which ones can be put to zero. 
Once we have the diagonal coefficients c^ 2k (r) we can 
calculate the off-diagonal ones with I < k with the help 
of equation (46). The coefficients with k < I < 2k are 
then clear for symmetry reasons. 

Using the computer algebra programme Maple V R7 we 
managed to calculate seven perturbative orders of the 
imaginary-time evolution amplitude which can be found 
at IE21. 



FIG. 2: This diagram shows which of the first-order coeffi- 
cients Cgujfr) have to be integrated (bold) and which ones can 
be obtained by employing symmetry and algebraic recursions 
(light). 



We have to expand the logarithm in order to obtain a 
perturbation expansion for the free energy F. For the 
first order we insert ( |l9| ) together with (||) into ( ff8| ) and 
evaluate the integral and for the second order we use, 
correspondingly, the data from Ref. ]l2|| . By taking the 
logarithm we get with (49) and with the expansion for 
the logarithm for the free energy to second order 



F( 2 )(/3) = I log 2 sinh ^ 



(50) 



AM 2 u. 



■ coth* 



hf3uj 



9 



2 h 3 



6AM 4 l 



^sinh 4 ^ 



36hf]u! cosh hj3u) + 60 sinh H0UJ + 21 sinh 2h(3u> 



sinh 



4 h/3ui 
2 



The higher orders are omitted for the sake of keeping the 
type face clear. 

With Maple we came as high as the fifth variational order 
which is two orders more than what has been obtained 
in previous work pM- 



VII. FREE ENERGY 

In this section we obtain perturbative results for the 
partition function by integrating the diagonal elements of 
our perturbative expression for the imaginary-time evo- 
lution amplitude from the previous sections: 



Z 



dx{x Kf3\x 0) . 



(48) 



From the partition function we then compute the free 
energy perturbatively: 



(49) 



VIII. DIAGRAMMATICAL CHECK 

It is possible to check the perturbative results for the 
free energy for all temperatures. Namely, we can expand 
Z in terms of harmonic expectations in a similar way as 
for the imaginary-time evolution amplitude in (|^). To 
that end we need the generating functional 



ZW)} 



+ OC 



dx(x hf3\x 0)uj [j] 



(51) 



which we get from (|7|)-(p0|). It is of the form 

ZW)] - z[o] 



(52) 



x exp 



1 

2fi2 
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where the harmonic partition function reads 

1 



Z[0] = 



2 sinh ; 



and 



G (p) (ri^ 2 ) 



h eosh ^ 



\Tl - T 2 U 



2Mu> 



sinh 



(53) 



(54) 



denotes the periodic Green's function of the harmonic 
oscillator. We now obtain the partition function Z of 
the anharmonic oscillator from the generating functional 
Z[j(r)] by differentiating with respect to the current j(r) 
while setting j(r) = afterwards: 



1 

exp< — - 



fit) 



drg 



Z[j(r)} 



3=0 



Thus we get 
Z = Z[0] 



(55) 



(56) 



hi) 



W 1 dTl 



hi) 



dr 2 



+ 72G^(t 1iTi )G^ 2 (n,T 2 )G^(T 2 ,T 2 ) 
+ 24G (p)4 (r 1 ,r 2 )] + ...} . 
In terms of Feynman diagrams this reads 



Z = Z[0] 




(57) 



,(58) 



= logZ[0] 



(59) 



where we have introduced the symbol 
1 
2 

Once we rewrite the partition function Z in the form of 
the cumulant expansion as on the right hand side of equa- 
tion ( p7| ) , the disconnected Feynman diagrams disappear 
|^| . Now we can easily take the logarithm. Following ( ff9| ) 
we obtain for the free energy 



F 




(60) 



The above Feynman diagrams are of course constructed 
with the help of the same rules as for the imaginary-time 
evolution amplitude (|l2|), (|l3|), and (|TJ), but instead of 
the Dirichlet's Green's function ([h]) we have to use the 
periodic Green's function (^4j). We now want to evaluate 
the four diagrams in (60) so that we get a second-order 
expression for the free energy for finite temperatures. Ac- 
cording to (|53| ) and ( |59|) we get for the zeroth-order con- 
tribution 



log 



2 sinh ^ 



whereas the first-order diagram becomes 



^ coth 2 ^ 



AM 2 



(61) 



(62) 



The integration in ( |62| ) is trivial, because G^ (r, r) docs 
not depend on r any more according to (|54|). For the 
second order the integrations become more sophisticated: 



ft^coth 2 -^ 
" 32M 4 w 5 sinh 2 5fii 
x (h(3uj + sinh hf3u>) . 

The other contribution to the second order yields 

256Af 4 w 5 sinh 4 ^ 
x (sinh 2h(3uj + 8 sinh hf3uj + 6h(3uj) . 



(63) 




(64) 
(65) 



So all in all we get for the free energy (|6C|) up to second 
order in the coupling constant g the result (^0[). It shows 
the correct low-temperature behaviour 



hm^ 2 >(/?) = ^ 



3gh 2 2lg 2 h 3 



4M 2 oj 2 8M 4 w 5 ' 



(66) 



which is the ground state energy and can be found for 
instance in |j| . 



IX. VARIATIONAL PERTURBATION THEORY 

Variational perturbation theory is a method that en- 
ables us to resum divergent Borel-type perturbation se- 
ries in such a way that they converge even for infinitely 
large values of the perturbative coupling ]|, pj| . To this 
end we add and subtract a trial harmonic oscillator with 
trial frequency f2 to our anharmonic oscillator (0): 



Mw 2 -n 2 



g- 



gx 



(67) 



Now we treat the second term as if it was of the order 
of the coupling constant g. The result is obtained most 
simply by substituting for the frequency to in the original 
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anharmonic oscillator potential 
ert's square-root trick 



id — » Sly/ 1 + gr , 



according to Klein- 



(68) 



where 



n 2 



gQ 2 



(69) 



These substitutions are not the most general ones. The 
square root is just a special case for the anharmonic os- 
cillator. 



We now apply the trick (68) to our first-order series rep- 
resentation for the free energy F found in ((5^). Substi- 
tuting for the frequency u> according to (pSj), expanding 
for fixed r up to the first order in g and re-substituting 
for r according to (|69|) we get 



1 



F(1)(M S l0S 2sinhIM 



am 2 q 2 



coth' 



+ 



1 coth 



(70) 



So the free energy ( |70| ) now depends on the trial fre- 
quency Q which is of no physical relevance. In order to 
get rid of it, we have to minimize its effect by employing 
the principle of least sensitivity |lg| . This principle sug- 
gests to search for local extrema of F((3, f2) with respect 
to ft: 



on 



(71) 



For the first order F^'(f3, f2) it turns out that there are 
several extrema for each (3. As we seek a curve Q,^ l \(3) 
that is as smooth as possible the choice is easy — we take 
the lowest branch for the others are not bounded (see 
Fig. ||). Moreover the other branches lead to diverging 
results. 

To second order, we proceed in a similar way and we 
find that there are no extrema at all for F^ ([3,Q). In 
accordance with the principle of least sensitivity we look 
for inflection points instead, i.e. we look for solutions to 
the equation 



d 2 F ( - 2 \f3,n) 
on 2 



= o. 



In general we try to solve the equation 







(72) 



(73) 



for the smallest possible n. Plugging QS N >(0) into 
F^((3,iY), we finally get back a re-summed expression 
for the physical quantity F{j3). The results for the first 
three orders are given in Fig. |4| In order to check our 
results we have to compare them to the numerically eval- 



uated free energy F^umiP) which is discussed in Sec. IXl 



p< 3 

S 2 

~G 

l 








P [l/2Ry] 



FIG. 3: Branch of the variational parameter fl^dJ) which 
we chose. The coupling strength is g = 1. Other branches 
not shown in this figure lead to highly diverging results. 
Throughout this paper all results are presented in natural 
units h = k B = 1 and, additionally, we have set M = w = 1. 



1,00 



^0,75 
Pi 

<N 

^0,50 
0,25 
0,00 



first order 




second and third order 
and numerical results 



P [l/2Ry] 

FIG. 4: Free energy of the anharmonic oscillator up to third 
order for intermediate coupling g = 1. The solid line repre- 
sents the numerical result FiS n(73h obtained by approximat- 
ing the partition function ( |7q ) with the help of the first ten 
energy eigenvalues. The other lines are variational perturba- 
tive results: The dashed line shows the first order, the dotted 
line shows the second order, and the dot-dashed line repre- 
sents the third order. Note that the second and third order 
are hardly distinguishable from the numerical results. Higher 
orders for a sp_ecial value of the inverse temperature /3 can be 
found in Fig. 



X. HIGHER ORDERS 

We now evaluate the convergence behaviour for the 
variational perturbative results for the free energy 
F^ N '(P) up to the fifth order. However, in order to re- 
duce the computational cost we restrict ourselves to a 
certain value of the inverse temperature f3. Results are 
shown in Fig. ^. For odd variational perturbation orders 
we optimized the free energy according to (f7l|), i.e. we 
determined f2 by setting the first derivative of F^ N>> ((3) 
with respect to f2 to zero. For even orders we had to go 
for inflection points, instead, so we had to solve equation 

It turns out that odd and even orders oscillate about 
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FIG. 5: The free energy of the anharmonic oscillator for in- 
termediate coupling g = 1 for f3 = 1 up to fifth variational 
perturbative order. The values converge exponentially to- 
wards the numerical value i ? iuJ n (l)=0.6571, which is shown 
as a straight line. The dimension of the free energy in natural 
units is 2Ry. 



an exponential best fit curve. For (3 = 1 the numerical 
result, J^um(l-O) = 0.6571, turns out to lie within the 
interval obtained by fitting the five perturbative orders 
of the free energy with origin, as shown in Fig. |5[ This 
interval is [0.657,0.659], and clearly, the variational per- 
turbative results converge exponentially. 
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Tl 







0.8037701932 


5 


14.203064494 


1 


2.7378891484 


6 


17.633934116 


2 


5.1792814619 


7 


21.236268598 


3 


7.9423804544 


8 


24.994705012 


4 


10.963538555 


9 


28.896941521 



TABLE I: The first ten energy eigenvalues E n of the anhar- 
monic oscillator for intermediate coupling g = 1. They were 
obtained by numerically integrating the Schrodinger equation 
with the initial condition that *(0) = 1, *'(0) = for even 
n, and ^(0) = and ^'(0) = 1 for odd n, and of course 
|^(ir)| < oo for large x. The energy eigenvalues are given in 
units of 2Ry. 
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XI. CHECKING OUR RESULTS 

The spectral representation of the partition function 
reads 

oo 

Z = £e-^», (74) 

n=0 

where the E n are the energy eigenvalues. Let us define 
the numerical approximants 

N 

^ = £e-^ (75) 

n=0 

and 

FW=-Ilog^W, (76) 

P 

respectively. One possibility to obtain numerical results 
for the eigenvalues E n is provided for by the so called 
"shooting method" . We integrate the Schrodinger equa- 
tion numerically for the potential (^|) and for a particular 
value of the coupling strength g. If the energy E which we 
plug into the program does not coincide with one of the 
energy eigenvalues E n , the solution to the Schrodinger 
equation explodes already for relatively small values of 
the coordinate x. If the energy eigenvalue is close to the 
exact answer, we have \^(x) \ < oo also for larger values of 
x. This method yields the unnormalized eigenfunctions 



FIG. 6: The numerial free energy i*num(/3), the first-order 
variational perturbative results for the free energy F^'(P), 
and the classical free energy F c i(/3) — — (log Z c i)/f3, from top 
to bottom. For small values of the inverse temperature j3 the 
classical calculations coincide with the other results. Lower 
temperatures, corresponding to higher values of j3, reveal dif- 
ferences between the classical approach ( |78| ) and quantum 
statistics. 



(the wave functions which still have to be normalized) 
and the energy eigenvalues to very high accuracy (see 
Tab. ||). Renormalization is necessary, for the computer 
algebra programme needs an initial value 'J'(O) which we 
set to one. Plugging the first ten numeric energy eigenval- 
ues into equation (|7|) and evaluating ((76|) up to N = 9, 

we obtain a function Fnum(P). It converges rapidly for 
low temperatures, corresponding to high values of j3. For 
high temperatures more terms should be taken into ac- 
count. 

Alternatively one can also use classical results as a high- 
temperature cross check: High temperatures correspond 
to classical statistical distributions such that we can eval- 
uate the classical partition function according to 



Z cl = / exp[-(3V(x)} , (77) 

J-oo A th 
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with the potential 
integral reduces to 



and A th = y/2wh 2 /Mk B T. This 



1 



'Mw 2 



x exp 



2A t h y 2g 

(3M 2 u 4 
32 5 



/4 



(3M 2 cj 4 
32g 



(78) 



where K 1 ^(z) is a modified Bessel function. The clas- 
sical partition function ([78]) can be evaluated for high 
temperatures which corresponds to small values of [3. 
Consequently, when we test our variational perturbative 
results, we compare them to the classical free energy for 
low values of (3, namely (3 < 1/4. And for high values 

of (3 we use the numerical approximation, -Fnum(/3), for 
comparison (see Fig. ||). 

In natural units h = k B = I a value of (3 = 1/4 corre- 
sponds to a physical temperature of T = 1.26 x 10 6 K. 



XII. CONCLUSION AND OUTLOOK 



turbation theory was found to be exponential. The fact 
that the principle of least sensitivity jl3| produces ex- 
trema for the odd variational orders and inflection points 
for even orders is reflected in the respective convergence 
behaviours: Odd and even orders can best be fitted sep- 
arately by exponentials as emphasized in Ref. fllf| . Thus 
we obtained intervals of convergence for certain values 
of the free energy which always turned out to contain 
the exact numerical result when taking into account the 
statistical errors associated with the boundaries of the 
intervals. For the free energy, the numerical results were 
obtained using its spectral representation reverting on 
the first ten energy eigenvalues obtained with the "shoot- 
ing method", sketched in Sec. XI. 



Finally, we note that our high-order perturbative results 
for the anharmonic imaginary-time evolution amplitude 
are useful for calculating other thermodynamic quanti- 
ties as the correlation function or the ground state wave 
function ||, |l5|, |l6| . Furthermore, it remains to compare 
these perturbative results with the semiclassical approx- 
imation |l7|| . 



The recursive techni que that has been developed 
throughout Sections IV -VI definitely out classes all di- 
agrammatical perturbative calculations. Using the con- 
ventional evaluation of Feynman diagrams, the partition 
function and the free energy have been evaluated up to 
third order jl4|, here we obtained the fifth-order result. 
For the free energy the convergence of variational per- 
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